#!/usr/bin/perl -w
# Aquanal, a perl script for the analysis of aqua output.
#
# This program is copyrighted by Kris Boulez (Kris.Boulez@rug.ac.be) and
# can be used under the license under which Perl is distributed.
#
# v 0.01 (11/10/95) :
# 	eerste versie
# v 0.02 (20/10/95)
#       Versie met ReadNrv package
# v 0.03 (1/12/95)
#       Veranderingen na test door Ton Rullmann
# v 0.04 (6/12/95)
#       verdere veranderingen na test van 0.03 door Ton Rullmann
# v 0.05 (22/4/96)
#	first version for release with aqua 0.40

BEGIN {
unshift(@INC,$ENV{'AQUASCRIPTSDIR'}); 
}
require 5.0;
use ReadNrv;
$ver = "0.05";
print "\n\nAquanal($ver): for the analysis of aqua output\n";
print "----------------------------------------\n\n";
print "If you want to use the default answer, which is printed between\n";
print "square brackets, just press enter for that question\n\n";
distrib (@ARGV);   # Which analysis to do 



sub anal_neighbour {  # $res1 , $res2, $mar (we zoeken +/-$mar)
shift;
use ReadNrv qw(@rv);
local ($r1,$r2,$t1,$t2,$r11,$r12,$r22,$r21,$t) = @_;
print "\nAnal_neighbour: analysis of which restraints there are between two\n";
print "residues (or ranges of residues). Ranges of residues are entered\n";
print "in the form (20-23). Overlapping ranges are not allowed.\n\n";
readnrv ($nrv,'g','rv');
print "The .nrv file contains $nres residues and $nrestr restraints\n\n";
$r1 = "" if (!defined ($r1));
$r2 = "" if (!defined ($r2));
do {
  print "First residue/range [$r1]: "; $t = <STDIN>; chomp($t);
  $r1=$t if ($t ne "");
  print "Second residue/range [$r2]: "; $t = <STDIN>;chomp ($t);
  $r2=$t if ($t ne "");
  print "Output to (s)creen/(f)ile [s]: "; $t = <STDIN>; chomp ($t);
  if ($t =~ /^f/) {
    print "  File name: "; $t = <STDIN>; chomp ($t);
    open (OUT, ">$t") || die "Can't open OUT $t: $!";
    select (OUT);
  }
  $r1 =~ s/\s*//g; $r2 =~ s/\s*//g;
  $t = 1;
  if ($r1 =~ /(\d*)-(\d*)/) {
    ($r11,$r12) = ($1,$2);
  }
  else {
    $r12 = $r1;
    $r11 = $r1;
  }
  if ($r2 =~ /(\d*)-(\d*)/) {
    ($r21,$r22) = ($1,$2);
  }
  else {
    $r21 = $r2;
    $r22 = $r2;
  }
  if ($r11 > $r12) {$t--;}
  if ($r21 > $r22) {$t--;}
  if (($r11 <= $r22)&&($r11>=$r21)) {$t--;}  
  if (($r21 <= $r12)&&($r21>=$r11)) {$t--;}
} until ($t > 0);
print "\n";
for (@rv) {
  ($t1,$t2) = (split)[2,6];
  if (inr($t1,$r11,$r12) && inr ($t2,$r11,$r12)) {next;}
  if (inr($t1,$r21,$r22) && inr ($t2,$r21,$r22)) {next;}
  if (inr($t1,$r11,$r12)||inr($t1,$r21,$r22)) {
    if (inr($t2,$r11,$r12)||inr($t2,$r21,$r22)) {
      local(@t,$t1,$t2);@t = (split)[0 .. 7];
      if ($t[0] eq "_") { $t[0] = ""};
      if ($t[4] eq "_") { $t[4] = ""};
      $t1 = $t[1].$t[0].":".$t[2].":".$t[3];
      $t2 = $t[5].$t[4].":".$t[6].":".$t[7];
      printf("%-12s- %-12s\n",$t1,$t2);
    }
  }
}
select (STDOUT); close OUT;
do {
  print "\n\n1. Run the same analysis again\n2. Main menu\n";
  print "\nq. Exit the program\n\n";
  print "Make a choice: "; $t = <STDIN>;chomp ($t);
} until ($t =~ /[12q]/);
if ($t eq "q") {exit;}
if ($t == 1) { anal_neighbour($nrv,$r1,$r2);}
if ($t == 2) { distrib ($nrv); }
} # anal_neighbour


sub inr {  # $t,$r,$m return 1 if  ($r-$m) <= $t <= ($r+$m) | anders 0
local ($t,$r1,$r2) = @_;
if (($t>=$r1) && ($t<=$r2)) {
  return (1);
}
else { return (0); }
} # inr

sub distrib {         # decides which function should be called
local ($t);
if ($#_ == -1) { $nrv = "";}
else { $nrv = shift; }
do {
  if (!-f $nrv) { $nrv = ""; }
  print "\nName .nrv file [$nrv]: ";
  $t = <STDIN>; chomp ($t);
  $nrv = $t if ($t ne "");
} while (!-f $nrv);
do {    # if there are args left we shouldnt do this
  print "\n\nThese are the different analysis\n\n";
  print "1. Which NOE's are there between certain residues\n";
  print "2. Highest restraints violations per structure\n";
  print "3. Average violation and range of violation values over all structures\n";
  print "\nq. Exit the program\n\n";
  print "Make a choice (1-3 q): ";
  $t=<STDIN>; chomp $t;
} until ($t =~ /[123q]/);
if ($t eq "q") {exit;}
elsif ($t==1) { anal_neighbour($nrv); }
elsif ($t==2) { max_viol_struc($nrv); }
elsif ($t==3) { max_viol ($nrv); } 
} # distrib

sub max_viol_struc {
shift;
local($nv,$k,$t,$r,$f,$file);
($nv,$k,$t) = @_;
$file = 0;
readnrv ($nrv,'g','rv');
use ReadNrv qw (@rv);
if (!defined ($nv)) {$nv = "";}
do {
  print "\nThere are $nrestr restraints and $nstruc structures\n\n";
  print "Nr. of violated restraints to print [$nv]: ";
  $t = <STDIN>; chomp ($t);
  $nv = $t if ($t ne "");
} until ($nv>0 && $nv<=$nrestr);
print "Output to (s)creen/(f)ile [s]: "; $t = <STDIN>; chomp ($t);
if ($t =~ /^f/) {
  print "  File name: "; $t = <STDIN>; chomp ($t);
  open (OUT, ">$t") || die "Can't open OUT $t: $!";
  $| = 1;select (OUT);
}
$t = select; if ($t eq "main::OUT") {$file = 1;}

for $k (1 .. $nstruc) {
  local (@t) = ();
  print "\nStructure $k\n-------------\n";
  sort_array ($k);
  ($t) = (split (/\s+/,$rv[0]))[$k+10];
  if ($t != 0.0) {
    print "\nLower bound violations\n----------------------\n";
    print STDOUT "." if ($file == 1);
    $f = 0;    # triggert of we al het aantal gevraagde restr geprint hebben
    for $t (0 .. $nrestr) {
      @t = ();
      $f = 1 if ($t == $nv);
      @t = (split (/\s+/,$rv[$t]))[0..7,$k+8];
      last if ($t[8] == 0); 
      last if ($f && $t[8] != $r);
      $r = $t[8];
      if ($t[0] ne "_") { $t[1] .= $t[0]};
      if ($t[4] ne "_") { $t[5] .= $t[4]};
      printf ("%6s %3d %-4s- %6s %3d %-4s  %.2f\n",$t[1],$t[2],$t[3],$t[5],$t[6],$t[7],$t[8]);    
    }
  }
  print "\nUpper bound violations\n----------------------\n";
  print STDOUT "." if ($file == 1);
  $f = 0;
  for ($t=$#rv;$t >= 0;$t--) {
    @t = ();
    $f =1 if ($t <= $#rv-$nv );
    @t = (split (/\s+/,$rv[$t]))[0..7,$k+10];
    last if ($t[8] == 0);
    last if ($f && $t[8] != $r);
    $r = $t[8];
    if ($t[0] ne "_") { $t[1] .= $t[0]};
    if ($t[4] ne "_") { $t[5] .= $t[4]};
    printf ("%6s %3d %-4s- %6s %3d %-4s  %.2f\n",$t[1],$t[2],$t[3],$t[5],$t[6],$t[7],$t[8]);
  }
} # for $k (1 .. $nstruc)
print STDOUT "\n" if (defined $file);
close OUT ;select (STDOUT); $| =0;
do {
  print "\n\n1. Run the same analysis again\n2. Main menu\n";
  print "\nq. Exit the program\n\n";
  print "Make a choice: "; $t = <STDIN>;chomp ($t);
} until ($t =~ /[12q]/);
if ($t eq "q") {exit;}
if ($t == 1) { max_viol_struc($nrv,$nv);}
if ($t == 2) { distrib ($nrv); }
} # max_viol_struc

sub sort_array {	# sorteren van een arrray (array , nr van de struct)

local ($k,$com,$l) = @_;
$k +=10;     # 9e kolom is eerste struc
$tmpfile = "/tmp/tosort.$$";
$tmps = "/tmp/sorted.$$";
open $tmpfile, ">$tmpfile";
for $l (0 .. $#rv) {
  print $tmpfile "$rv[$l]\n";
}
close ($tmpfile);
$com = "sort -o ". $tmps . " +".$k."n ". $tmpfile;
system ("$com"); # the actual sort goes via the system sort.
open TMP, "$tmps";
@f = ();
while (<TMP>) {
  chomp;
  push @f,$_;
}
@rv = @f;
unlink $tmpfile;
unlink $tmps;
}


sub max_viol {  # max_viol kijkt of er in 1 van de $nstruc een viol optreedt
                # die groter is dan max_viol. In dat geval wordt de gem. viol,
                # stand. dev, max- en minimum viol bepaalt

local ($k,$t,@t,$max,$min,$mir,$mar,@temp,$tot,$avg,$std);
shift;
local ($max_v) = @_;
use ReadNrv qw(@rv);
readnrv ($nrv,'g','rv');
print "\nThere are $nstruc structures and $nrestr restraints\n\n";
print "Each restraint which is violated in one of the $nstruc structures by\n";
print "more than the maximum violation will be listed.\n";
$max_v = "" if (!defined $max_v);
print "\nMaximum violation [$max_v]: ";
$t = <STDIN>; chomp ($t);
$max_v = $t if ($t ne "");
print "Output to (s)creen/(f)ile [s]: "; $t = <STDIN>; chomp ($t);
if ($t =~ /^f/) {
  print "  File name: "; $t = <STDIN>; chomp ($t);
  open (OUT, ">$t") || die "Can't open OUT $t: $!";
  select (OUT);
}
print "\n  Atom 1      Atom 2         Average viol.   Individual viol.     Restraint\n";
print "                             Ave.    Std.  Low (str)  High (str)  Low  Upper\n";  
print "-" x 77; print "\n";
for $k (0 .. $#rv){
  @t = split (/\s+/,$rv[$k]);
P: for $t (11 .. $nstruc+10) {     # eerst checken of we iets
     if ($t[$t] > $max_v) {          # moeten doen
        $max=0;$min=1000;$tot=0;$std=0;         # beginwaarden
        for $t (11 .. $nstruc+10) {   # opnieuw beginnen tellen
          $tot += $t[$t]; $temp[$t-11]=$t[$t];
          if ($t[$t] > $max) {$max=$t[$t];$mar=$t-10;}
          if ($t[$t] < $min) {$min=$t[$t];$mir=$t-10;}
        }
        $avg = $tot/$nstruc;
        for $t (0 .. $nstruc-1) {
          $std += ($avg-$temp[$t])**2; }
        if ( $nstruc == 1 ) {
            $std = sqrt ($std/$nstruc);
        } else {
            $std = sqrt ($std/($nstruc-1));
        }
	if ($t[0] ne "_") { $t[1] .= $t[0]};
	if ($t[4] ne "_") { $t[5] .= $t[4]};
      printf("%5s %2d %-4s %-6s %2d %-4s %.2f +/-%.2f [%.2f(%2d) - %.2f(%2d)]  (%.2f-%.2f)\n",
$t[1],$t[2],$t[3],$t[5],$t[6],$t[7],$avg,$std,$min,$mir,$max,$mar,$t[9],$t[10]);
        last P;
        }
     }
  }
close OUT ;select (STDOUT);
do {
  print "\n\n1. Run the same analysis again\n2. Main menu\n";
  print "\nq. Exit the program\n\n";
  print "Make a choice: "; $t = <STDIN>;chomp ($t);
} until ($t =~ /[12q]/);
if ($t eq "q") {exit;}
if ($t == 1) { max_viol($nrv,$max_v);}
if ($t == 2) { distrib ($nrv); }
} # max_viol













